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Abstract 

We present additions to the martini event generator for examining heavy quarks and quarkonia 
in heavy-ion collisions. All stages of a heavy-ion collision affect the observables associated with 
heavy quarks: the initial phase space of the heavy quarks are sampled with pythia8.1, the heavy 
quarks are evolved using Langevin dynamics and a 3+1-dimensional hydrodynamical description 
of the heavy-ion collision, and are fragmented and hadronized using a modified version of the 
color evaporation model that takes into account non-trivial evolution in position space, as well 
as the possibility of recombinant quarkonium production in heavy-ion collisions. We use this to 
re-examine the production of quarkonium at RHIC, and anticipating vertex detection we predict 
yields of B c mesons at RHIC and the LHC. 
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I. INTRODUCTION 



Heavy quarks and quarkonia are hard probes of heavy-ion collisions; they are produced 
in hard processes in the initial moments of the collision and exist at all stages of the event. 
As a result, any calculation related to heavy quarks must consider nearly every aspect of 
heavy-ion physics: shadowing and anti-shadowing of the parton distribution functions, the 
perturbative cross sections, the finite-temperature evolution of heavy flavor, and finally, in- 
medium hadronization. Any calculation which neglects any of these steps and compares 
with an observable must quantitatively justify why the step was ignored. 

As a result, the heavy quark and quarkonium observables in heavy-ion collisions now have 
a history beginning with qualitative understanding followed by scrutiny with simulation. As 
an example, consider the changes in J/tp yields caused by finite-temperature effects. These 
changes were first described as being caused by the non-Abelian analog of the photo-electric 
effect [H [2] , and were estimated to be small. A large J / ip suppression was also predicted be- 
cause of Debye screening of the QQ potential, which had been seen in lattice calculations [3] ; 
the result suggested so large an effect that J/ip suppression was hoped to be the "smoking 
gun" for formation of quark-gluon plasma in heavy-ion collisions. An "anomalous" suppres- 
sion was indeed observed by the collaborations at the SPS [HE]. However, the results from 
RHIC indicate a more complicated suppression pattern than what can be described with the 
original results for large suppression the suppression at RHIC was small in comparison 
to expectations based on the results from SPS. Since then, heavy quark observables have 
been examined more carefully than had been done in the first semi-analytic predictions: 
the flow of heavy flavor has been estimated using hard thermal loop effective theory and a 
hydrodynamical description of heavy-ion collisions [7] ; the rates for destruction and regen- 
eration of J I ip particles at finite-temperature have been calculated using finite-temperature 
heavy quark effective theory and potentials from lattice QCD, and these rates were applied 
to the results at RHIC and the LHC [SHTU]: and finally, J /if) suppression was described with 
heavy quark diffusion and this description was integrated with a 2+1-dimensional hydrody- 
namical simulation for estimating changes in J / ip yields at RHIC [HI [12] . Most recently, the 
rates for thermal quarkonium suppression were calculated away from thermal equilibrium, 
for viscous plasmas, and applied to the suppression of bottomonium states at the LHC [13J. 

Event generation is a promising approach to keeping all of the steps under control while 
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adding new physical considerations. In this paper, we describe modifications made to the 
martini event generator [T3], which has been used successfully for examining jets and pho- 
tons at RHIC and the LHC [ToT - ITT] . We cover the physics of heavy quarks "chronologically" : 



in Section II A production, then in-medium evolution in Section II B, and finally hadroniza- 



tion in Section |II C| The new physical considerations of this paper are mostly in the section 
on hadronization, and their implementation in MARTINI represents an important part of 
what Monte Carlo can do that is difficult analytically without many simplifying assump- 



tions. We compute the results for observables at RHIC and the LHC in Section III and 
make a prediction of the yields of B c mesons, which are produced in heavy-ion collisions 
almost entirely recombinantly. Finally in the conclusions, areas where future work is needed 
are pointed out. 



II. THE STEPS OF HEAVY QUARK EVENT GENERATION 
A. Initial production of heavy quarks 

Because of the masses of charm and bottom quarks, which are large compared with Aqcd 
and with the temperatures accessible in heavy-ion collisions, the processes producing these 
quarks have a large energy scale and can be described perturbatively. Heavy quarks are like 
high-energy jets in that their initial production in a heavy-ion collision is some elaboration 
of the production in proton-proton collisions. Also, the thermal production of charm and 
bottom quarks in heavy-ion collisions is likely small: t-channel gluon-gluon scattering which 
contributes most to the thermal heavy quark production is suppressed by two factors of the 
initial Bose-Einstein distributions of these gluons. For this reason, we consider only heavy 
quarks produced initially in the nuclear collision. 

However, even perturbative production of heavy quarks is complicated by nuclear ef- 
fects in heavy-ion collisions. First of all, the parton distribution functions in a heavy- 
ion collision are modified because of isospin. Every nucleon's parton distribution can be 
described with the weighted average of the neutron and proton distribution functions, 
fz(x, Q) = ^f p (x, Q) + ^j^-f n (x, Q), to take this effect into account. An effect with greater 
significance for heavy quark production is the shadowing of the parton distribution func- 
tions by nucleons at the same point in the transverse plane, and the effect of anti-shadowing 
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required by momentum conservation. These modifications have a satisfying description in 
color glass condensate effective theory, and if the saturation scale Q s > 2rriH the color glass 
condensate would be the best tool for describing heavy quark production at RHIC and the 
LHC [IH]. Gluon saturation likely does cause some of the anomalous nuclear modification 
of J ftp production, however it cannot explain all of the observed modification. Another 
approach is more phenomenological and practical for event generation: parametrizing the 
parton distribution functions for nucleons in bound states, and simply evolving these ac- 
cording to leading-twist DGLAP evolution [2U] . 

Another complication in the initial production arises at high-p^ relative to the heavy 
quark's mass. In the factorized approach to heavy quark production, the perturbative 
production at next-to-leading order contains terms depending on logfj^r/mg), which need 
to be resummed when px ^> ttiq. Working with "massless" fragmentation into hadrons can 
give the correct results at high p? while now neglecting terms proportional to \og(mQ / pt) ■ 
Cacciari et al. effectively merged the results at fixed order in perturbation theory (accurate 
at low p T ) and next-to-leading log order (accurate at high p T ) [2Tj . 

We use PYTHIA8. 1 to sample the six-dimensional momentum distributions of QQ systems. 
This works only at leading-order in perturbation theory, as almost all state-of-the-art event 
generators do not include higher orders; including all of the elements in pythia at higher 
orders is problematic because of increases in the complexity of phase spaces. However, some 
part of the higher-order behavior is considered by setting the total cross section of charm 
by hand (we use 0.6 [ib for proton-proton collisions at RHIC and 3.0 fib for the LHC). The 
nuclear shadowing is taken into account using EKS98 nuclear parton distribution functions, 
and isospin effects are taken into account by sampling from p+p, n+p, p+n, and n+n 
collisions. 

Finally, the spatial dimensions of the phase-space distribution for the QQ states must be 
sampled. We use the Glauber model to determine the initial centers of mass of the pairs 
in the transverse plane. An important additional consideration is required by the non-zero 
thermalization time To in heavy-ion collisions. We evolve the heavy quarks during this time 
according to the equations of motion determined by the zero-temperature Cornell potential: 

dx l • . dp 1 ^Cornell /, x 

~dt =p/E > = — (1) 

These equations are not manifestly Lorentz-invariant due to the frame-dependent Vcomdb 
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however, it is defined in the center-of-mass frame of the pair for reasons in Section II C 



making the equations of motion Lorentz-invariant while neglecting radiative effects. The 
equations of motion in the following section generalize these for temperature-dependent 
potentials and diffusive effects. The leapfrog method is used to achieve sufficient numerical 
accuracy: simpler discretizations of the equations of motion were found to have unacceptably 
large rounding errors at forward rapidities. 



B. Heavy quarks at finite temperature: Langevin dynamics 

Radiative "splittings" of high-energy partons, which are a significant source of in-medium 
evolution for high-pr massless partons, are suppressed for heavy quarks because of the 
"dead-cone" effect, already known from jet physics [22]. The evolution of heavy quarks at 
finite temperature is approximated as being caused by elastic collisions in the hard thermal 
loop (HTL) approximation of QCD. At the temperatures accessible in heavy-ion collisions, 
thermal heavy quarks exist in the "diffusive limit" of dynamics. To understand qualitatively 
why this is so, consider a heavy quark with mass M, with M 3> T. In the hard thermal 
limit of QCD, the heavy quark suffers kicks leading to momentum transfers on the order of 
gT. The typical momentum transfer is therefore much smaller than the heavy quark mass, 
meaning that it takes many kicks to change significantly the momentum of the heavy quarks. 
The dynamics of heavy quarks undergoing elastic processes can thus be approximated as 
diffusive, and therefore differential in time. 

Many general properties can be calculated directly from the Langevin equation. In the 
rest frame of the medium, 

= [(<% -PiPj)tor(p) +PiPj K L(p)} , 

(2) 

where 77 is a drag coefficient, £j the stochastic force on the heavy quark, and Kl {kt) the 
longitudinal (transverse) momentum diffusion coefficient. Often we will approximate k = 
k-t = ^Li which is true non-relativistically. These calculations, along with the requirement 
that an ensemble described by this equation thermalizes, leads to the Einstein relation 
77 = k/2MT. However to simulate individual events, being able to calculate ensemble 
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averages only is inadequate; this stochastic differential equation needs to be discretized. 
The central limit theorem provides the means to do this: over a sufficiently long timestep 
5t, the sum of many random kicks is described by a Gaussian with width a/ n5t. In the end, 
we use Monte Carlo also in solving the in-medium evolution equation. 

In the diffusive approximation, one physical parameter needs to be determined: Dhq, the 
spatial diffusion coefficient for the heavy quark. It can be related to k and 77 by examining 
the expectation value of ((x(t) — x(0)) 2 } = QDt with the Eq. Q, and can be calculated 
in the non-relativistic limit as caused by a t-channel gluon exchange in the hard thermal 
approximation of QCD [7J, which has been done at next-to-leading order [23] . It has also 
been calculated in supersymmetric large-iV c QCD at strong coupling using gauge-gravity 
duality |24j, as well as in gravity models of the dual for QCD [25] . The gravity calculations, 
as well as the poor convergence of HTL effective theory, point to a small diffusion coefficient, 
~ 1/2ttT. In the simulations that follow we set 27tTDhq = 3.0 and do not consider here the 
effects of tuning this parameter, which has large effects on open heavy flavor but smaller 
effects on the total quarkonium yields. 

The dead-cone effect causes brehmsstrahlung to be suppressed in the cone within an angle 
Pt/E around the direction of a heavy quark's momentum. We would expect this cone to 
shrink at high momenta, causing the dead-cone effect not to be significant and the effect 
of brehmsstrahlung to be no longer negligible, and indeed, the momenta accessible at the 
LHC are well into this range. For this reason, a future paper will determine the effect of 
brehmsstrahlung on heavy quarks, as well as describe the heavy quarks at lower momenta 
without the diffusive approximation. 

There is one final consideration for the dynamics of heavy quarks that we need to make: 
the self-interaction of QQ color singlets at finite temperature is not negligible, as can be 
seen with reliable lattice calculations of the correlation between two Wilson lines at finite 
temperature [25]. The effect of temperature on the potential cannot be neglected; in fact, 
the earliest prediction of significant J/ip suppression was based entirely on changes in this 
potential with temperature. 

This raises another physical question: what thermodynamic potential of QQ color singlets 
is appropriate to use, the internal energy used for adiabatic processes, or the free energy used 
for isothermal processes? This question has been investigated by comparing potential model 
results for quarkonium correlation functions with the results from lattice QCD [27J EE] • It 
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remains difficult to distinguish which potential is best for quarkonium at finite temperature. 
We continue to use the internal energy as a function of temperature, as was done in [TTl IT2] . 
based on the reasoning in [2H] comparing the timescales of binary bound states, related to the 
binding energies of these states, to the relaxation times of the medium at finite temperature. 

C. Heavy quark hadronization 

Heavy quarks hadronize into mesons with either one or two heavy quarks. Because 
heavy quarks are produced almost entirely in flavor-conserving processes, heavy quarks are 
produced with an anti-quark of the same flavor and both of these possibilities for hadronic 
final states need to be considered, particularly at center-of-mass energies near thresholds for 
quarkonium; the production of a given hadronic state should never be considered singly. 

First, consider open heavy flavor: a single heavy quark Q with momentum P fragments 
into a hadron Qq with momentum zP, leading to an energy difference AE between the 
initial and final state. General quantum-mechanical arguments lead one to estimate the 
probability for a given transition to be oc 1/(AE) 2 , leading to the Peterson fragmentation 
formula [30J: 



In MARTINI, whenever heavy quarks are determined to lead to open heavy flavor mesons, 
this function is sampled to determine the momentum fraction of the resulting meson. At 
this point in time, for simplicity, the various open heavy flavor states are not distinguished. 

The previous argument assumed that the heavy quark was sufficiently energetic (and 
slightly off-shell) so that a final open flavor state was kinematically accessible. Considering 
the possibility that this state is not accessible allows one to estimate the production of 
hidden heavy flavor states. Given a heavy quark and anti-quark, with momenta p\ and P2, 
respectively, whose spatial separation in the lab frame is given by r^, we define the modified 
invariant mass of the pair to be 



where V(r) may be reasonably approximated at chemical freezeout with the Cornell poten- 
tial. This invariant mass determines what final states are kinematically accessible to the 




(3) 




(4) 
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pair if it hadronizes adiabatically; if M12 < 2m^+ , open heavy flavor would be kinematically 
inaccessible to the pair and the final state must be quarkonium. 

In this way, our method of hadronization resembles the color evaporation model [HTl 
[32] . which examines simply the invariant mass of the pair to determine whether or not 
quarkonium is formed. The color evaporation model is quite successful in determining the 
quarkonium yields in experiments using a range of targets and projectiles, however, it makes 
no predictions concerning the important observable of J/if> polarization. The color singlet 
model [33] and applications of NRQCD [3UI3S] make accurate predictions of the J /if) yields 
only with some tuning but after being calculated to higher orders in perturbation theory are 
yielding better agreement while describing polarization. 

For the purposes of this paper, the distinctions between the mesons are often ignored. 
Specifically, Peterson fragmentation leads to the production of D mesons, without D + , D~ , 
and D* mesons distinguished, and while the color evaporation model can use the invariant 
masses of the pairs to separate the quarkonium into various excited states, there is no 
distinction made between i] c and the three polarizations of J /if). These states can have 
similar masses but different quantum numbers leading to different decay channels. For 
this reason, the final yields calculated with martini need to be multiplied by well-defined 
correction factors before being compared with data, for example the yields of ground-state 
quarkonium, when being compared with J /if) yields, needs to be multiplied by a factor of 
3/4 so that the r\ c particles are not included in the yields. For some processes, for example 
e + e _ annihilation, such a factor would be inappropriate (the r] c has the wrong quantum 
numbers for this process and is not produced at all), however in heavy-ion collisions these 
quantum numbers are best estimated as entirely randomized, leading to each state's yields 
being proportional to its degeneracy. 

Heavy-ion collisions lead to the possibility of recombinant quarkonium production. Be- 
cause of the masses of heavy quarks, there is at leading order only one hard process leading 
to heavy quark production per proton-proton collision; the production of multiple heavy 
quarks is subleading in a s . However in a heavy-ion collision, there can be many such hard 
processes: in a typical central gold-gold collision at RHIC ~20 heavy quark pairs are pro- 
duced and at the LHC ~100 pairs. This forces us to consider the possibility that a heavy 
quark from one hard process and an anti-quark from a different hard process might have an 
invariant mass below the threshold of 2rri£)+ and may form quarkonium. 
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To determine recombinant production, we do not assume that the heavy quarks hadronize 
adiabatically; instead the bulk of the collision is assumed to act like a heat bath at tem- 
perature T ch , the chemical freezeout temperature for quarkonium. The possible final states 
for the heavy quarks are assumed to be open heavy flavor mesons and quarkonium (baryons 
with open heavy flavor, while interesting, are subleading and not considered here). Each 
possible "pairing" of the 2N heavy quarks and anti-quarks is determined by some a, an 
element of the permutation group SV- Finally, we associate with each of these pairings an 
energy: 

N 

E a = Y,M ia{i) , (5) 

i 

where Mj a ^ is determined as in Eq. (|4j). Now that both an energy for each of these pairings 
and a temperature are defined, we can simply weight each pairing by its Boltzmann factor 
to determine the yields for recombinant production. 

The previous argument applies well for an infinite volume cooling uniformly, however in 
heavy-ion collisions, this limit is not appropriate for estimating recombinant quarkonium 
production. Heavy quarks reach their chemical freezeout surface at different times, and it is 
appropriate for a heavy quark Q to hadronize based on the phase space of heavy quarks at 
the time Q hadronizes. For this reason, the hadronization of heavy quarks is determined in 
MARTINI according to this three step algorithm: 

• When T < T c h for a heavy quark Q, the modified color evaporation model is used to 
determine whether or not it forms quarkonium with the other heavy quark created in 
the perturbative process that created Q. If quarkonium is formed, then the hadroniza- 
tion stops here. This channel for quarkonium production is possible for heavy quark 
pairs produced in regions of the transverse plane at temperatures both above and be- 
low the deconfinement phase transition, and takes into account the heavy quark pairs 
that are not significantly affected by in-medium dynamics. 

• If quarkonium has not formed, then recombinant production is considered. All op- 
positely charged heavy quarks are checked with the color evaporation model for re- 
combinant production with Q and all such heavy quarks are added to the list. If 
the set is non-empty, then one of these heavy quarks is sampled and a quarkonium 
state is formed with Q. In this situation the Boltzmann factor for open heavy flavor 
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production is so small that we neglect it. This channel for production is only available 
to heavy quarks that exist in QGP. This channel is also not available for the directly 
produced quarkonium (quarkonium formed in the previous hadronization step) as it 
is approximated as kinematically unavailable to the heavy quarks in these states. 

• If neither surviving or recombinant quarkonium has formed, the Peterson fragmenta- 
tion function for Q is sampled and an open heavy flavor meson is formed. 

Notice that we use two thermodynamic limits for heavy quark hadronization: adiabatic 
for Peterson fragmentation and color evaporation for direct production, and isothermal for 
recombinant production. While the use of these different limits is well-motivated physically, 
more work in this area can and should be done, especially since all hadronic observables are 
sensitive to this final step. 

III. RESULTS FOR RHIC AND THE LHC 

We now apply martini to observables at RHIC and the LHC. To maximize the un- 
derstanding of these effects, we will not consider feed-down of excited quarkonium or B c 
mesons into ground states of quarkonium. Vertex detection will allow semi-leptonic weak 
decays of B c and other b-hadrons to be observed relatively easily, and these hadrons to 
be identified, while the use of tracking should allow di-particle invariant mass distributions 
besides the dielectron invariant mass spectrum to be used to find excited quarkonium; for 
example, the 7 — J/ip invariant mass spectrum can be used to determine the yields of \ c in 
heavy-ion collisions. Indeed, one of our main points in this paper is that this analysis should 
be undertaken, considering the high statistics for charm yields at the LHC. 

A. Prompt and non-prompt J ftp suppression 

Fig. [T] shows the results of martini following from all of the above considerations, 
applied to gold-gold collisions at y/s NN = 200 GeV. The only cut applied to the final- 
state particles is that the absolute value of its pseudorapidity \r]\ < 0.35. Raa denotes the 
"anomalous" nuclear modification factor: a significant suppression is caused by inelastic 
scattering of quarkonium with the nucleons in the original gold nucleus, and this is not 
taken into account. 
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FIG. 1: The anomalous suppression of charmonium in Au+Au collisions at RHIC energies, with 
cuts in pseudorapidity appropriate for the PHENIX detector (\rj\ < 0.35). Note the sequential 
behavior of the suppression, and how event generation leads to larger recombinant production 
than is predicted by semi-analytic results. 

The pattern of suppression in these plots is consistent with sequential quarkonium disso- 
ciation: excited quarkonium is destroyed more efficiently than is quarkonium in the ground 
state, due to the ground state's large binding energy relative to the temperature [36] • How- 
ever, this suppression pattern appears here without enforcing the sequential dissociation 
above predetermined dissociation temperatures for the quarkonium states by hand; this 
pattern occurs simply because less energy has to be transferred from the medium to a QQ 
system with a small binding energy to excite it above the 2m d threshold than needs to be 
transferred to a more more tightly-bound QQ system. Also noteworthy is the size of the 
recombinant contribution to production. It is larger than was estimated in [12J. Simple es- 
timates of recombinant production fail to take into account the strong correlation of charm 
quarks in position, because of their scaling with number of collisions in the transverse plane 
which peaks sharply at the center. Also, martini, by simulating the full event, correctly 
fluctuates the number of charm quarks per event and considers events with more charm 
quarks than average, taking into account the large recombinant production possible in those 
events. 

Fig. [2] shows the same algorithms applied to lead-lead collisions at yfs NN = 2.76 GeV. 
Recombinant production of excited charmonium leads to an increase in Raa with A^ part , due 
to large charm densities at small rapidity. Note especially how including recombinant pro- 
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FIG. 2: The anomalous suppression of quarkonium for Pb+Pb collisions at LHC energies, with 
cuts appropriate for the CMS detector. 

duction of excited states leads to a larger Raa than that of the J/tp yields. The recombinant 
production explains why the sequential suppression pattern was not observed by CMS. We 
predict that if the nuclear "absorption" cross-section for charmonium is taken into account 
properly as in [37], the CMS results could use the proper baseline for the initial quarkonium 
production and would show this increase in non-prompt J ftp Raa- Such a characterization 
of very many cold nuclear matter effects can best be understood with accompanying pA 
data. 

Using martini, heavy quark flow can be re-examined. Fig. [3] shows V2(j>t) for charm 
both before and after hadronization, and for two kinetic freezeout temperatures, Tf Q = 
190 MeV and Tf Q = 135 MeV. The v 2 result for open charm with the high kinetic freeze-out 
temperature is smaller than the result in [7J for 2ttTDuq = 3 by about a factor of one half, for 
a simple reason: the flow of heavy quarks develops slower than the flow of the bulk, so that 
quicker freeze-out leads to smaller azimuthal anisotropics. The smaller kinetic freeze-out 
temperature for open charm leads to larger azimuthal anisotropy, and therefore, sequential 
kinetic freeze-out can explain the different results for the flow of open and hidden charm. A 
simple argument related to the binding energy of mesons, as well as the results from gauge- 
gravity duality [38J, suggest that the kinetic freezeout temperature for quarkonium should 
be larger than the freezeout temperature for open heavy quark mesons. In this way, our 
plot shows how a small ^(pt) for quarkonium, relative to the t>2 for open heavy flavor, can 
be explained. Fig. [4] shows our results for heavy quarkonium flow at the LHC, compared 
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FIG. 3: V2(pt) for charm mesons in a Au+Au collision with y / s/v r /v = 200 GeV, when 6 = 5 fm. 
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FIG. 4: V2(pt) for charm mesons in a Pb+Pb collision with ^/snn = 2.76 TeV, when 6 = 5 fm. 

with the flow of single charm mesons at the same kinetic freezeout temperature. Finally, 
note the important effect of fragmentation on the ^(pt) results for open charm mesons. 

B. B c mesons 

Schroedter et al. predicted B c meson production in heavy-ion collisions, estimating the 
production to be one such meson every 20 collisions [39|. Detecting these states would be 
difficult with the results from the original incarnations of the detectors at RHIC: the vertices 
of their decays were not detected, instead the momenta of the relatively stable results of 
their weak and electromagnetic decays are measured, analyzed, and binned according to pair 
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FIG. 5: Predicted B c meson yields at both RHIC and LHC energies, in the pseudorapidity ranges 
of the PHENIX and CMS detectors, respectively, as a function of centrality. 

invariant masses to determine their yields. This makes the progenitors of two-body decays 
(such as the J/ip) difficult enough to observe, and makes the observation of B c mesons by 
their three-body decay challenging. 

This situation should change with vertex detection at RHIC as well as at the LHC. The 
capability of CMS in this regard had already been demonstrated in distinguishing x from 
J/i/j mesons in proton-proton collisions. To date, the results for quarkonium yields in heavy- 
ion collisions were made using the dielectron spectra from these collisions; this leads to all 
results suffering from the feeddown caused by excited states and makes the determination 
of the yields of particles which do not have simple decays into dileptons challenging at best. 

Fig. [5] shows martini's prediction for B c meson yields, at both RHIC and LHC energies, 
for the pseudorapidity ranges detectable by the PHENIX and CMS detectors, respectively. 
It is worth noting that according to these predictions, heavy-ion experiments can be the 
main means for studying these rare states. Also, as these states are produced entirely 
recombinantly, they are a unique probe of finite-temperature hadronization. 

IV. CONCLUSIONS 

Using full event generation of heavy-ion events, we have examined and re-examined heavy 
quark observables at RHIC and the LHC. We have worked completely in the diffusive limit 
of the heavy quark dynamics: the quarks were described with Langevin equations of motion 
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with a potential term, derived from lattice QCD simulations. Hadronization involved an 
elaboration of the usual color evaporation model. Full event generation allowed for reliable 
estimation of recombinant quarkonium production, including the production of B c mesons. 
We found that this recombinant production is larger than previous estimates. 

The first improvement that needs to be made here is the elimination of the "diffusive ap- 
proximation" for heavy quarks. In order to do so, we can calculate the differential HTL rates 
dT(p)/d 3 q and sample these quite easily. However, the radiative rates are also important for 
heavy quarks at the LHC, where transverse momenta are accessible for which the dead-cone 
effect is less significant. Work in progress is eliminating the diffusive approximation from 
these simulations. 
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